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We introduce a Generalized Loop Move (GLM) update for Monte Carlo simulations of frustrated 
Ising models on two-dimensional lattices with bond-sharing plaquettes. The GLM updates are 
designed to enhance Monte Carlo sampling efficiency when the system's low-energy states consist 
of an extensive number of degenerate or near-degenerate spin configurations, separated by large 
energy barriers to single spin flips. Through implementation on several frustrated Ising models, 
we demonstrate the effectiveness of the GLM updates in cases where both degenerate and near- 
degenerate sets of configurations are favored at low temperatures. The GLM update's potential 
to be straightforwardly extended to different lattices and spin interactions allow it to be readily 
adopted on many other frustrated Ising models of physical relevance. 
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I. INTRODUCTION 

Monte Carlo (MC) simulations are among the most 
ubiquitous computational tools used in statistical physics 
and material science. Modern MC methods, evolved from 
the original Metropolis method of the 1950s [lj, have be- 
come increasingly sophisticated with the advent of cluster 
moves, histogram reweighting methods, parallel temper- 
ing, and other technical advances [2(- This sophistica- 
tion, coupled with the continuing increase in available 
computing power, allows MC methods to simulate classi- 
cal and quantum statistical mechanical systems of a level 
of complexity unimaginable even a decade ago. Indeed, 
the remarkable growth in the size of systems accessible to 
MC simulations is owed as much to the advancement of 
algorithm technology as it is to the increase in available 
raw CPU power through Moore's law Q. 

Lattice magnetic systems offer some of the greatest 
challenges to MC practitioners. In fact, it was realized 
early on that MC simulations of the ferromagnetic Ising 
model employing simple local updates (single spin flip 
(SSF) algorithms) suffer severe critical slowing down - a 
rapid increase in autocorrelation times - near a second- 
order phase transition. This stimulated work on non- 
local (cluster or collective-mode) algorithms, such as the 
well-known Swcndscn-Wang [4[ and Wolff [5J algorithms. 
A different situation that also requires non-local updates 
is known to occur in broad classes of frustrated magnetic 
Ising models. Such models, typically identified by pre- 
dominantly antifcrromagnetic (AFM) or random inter- 
actions, can have disordered groundstates which consist 
of an extensive number of equal-energy (degenerate) spin 
configurations - such is famously the case for the trian- 
gular lattice Ising model with AFM interactions jy, 0j- 
In such models, simple SSF updates have the tendency 
to bring the configuration out of the degenerate manifold 
of groundstates, hence costing an energy proportional to 
the Ising interaction. At sufficiently low temperatures, 
the Metropolis algorithm will reject such moves, inhibit- 
ing ergodicity. which is the ability of the MC simulation 



to explore the entire degenerate manifold of states in rea- 
sonable (i.e. non-exponential) computing time. 

This issue was recently brought to the forefront in the 
broad class of ice or vertex models [a] . Of largest interest 
are the so-called spin ice models [9j - Ising models on a 
frustrated pyrochlorc lattice, which are realized experi- 
mentally in some rare-earth titanate compounds. Ideal 
spin ice Ising models promote a disordered degenerate 
groundstate. However, it was found that weak pertur- 
bations, such as occur from long-range dipolar interac- 
tions, can energetically favor one or several configurations 
[101 ] . Exploration of different configurations contributing 
to the degenerate groundstate is thus a crucial task for 
MC simulations, and can only be achieved through loop 
updates [II], [Ijl • The dynamics of loops in spin ice and 
related models are also intimately tied to important phys- 
ical phenomena such as the Kasteleyn transition [1 31 ] , and 
the concept of fluctuations between topologically ordered 
ground-state sectors [l4[. In addition, loop dynamics in 
MC simulations have become a crucial part of our un- 
derstanding of aspects of experimental dynamics in real 
spin ice materials [l2J . 

Since the identification of loop algorithms is topical 
not only as a simulation technique but as a probe into 
the physics of materials, it is remarkable that generalized 
loop algorithms do not exist for broad classes of different 
frustrated models. Rather, each system typically needs 
to have an algorithm designed based on the particular 
constraints of its groundstate manifold. Such is the case 
for vertex models @, spin ice [Il],[l2| ( an d related corner- 
sharing triangle models such as the Ising kagomc AFM), 
and recently the fully-frustrated honeycomb Ising model 
[l5| . In the present work, we address this shortcoming by 
introducing a generalized loop algorithm for a large class 
of frustrated Ising models with bond-sharing plaquettes. 

The paper is organized as follows. In Section [Hi we 
outline a Generalized Loop Move (GLM) update in de- 
tail for a class of two-dimensional frustrated Ising mod- 
els, illustrating the specific examples of the triangular 
lattice Ising AFM, and the fully-frustrated square and 



honeycomb lattice Ising models. We prove rigorously 
that the algorithm obeys detailed balance (Appendix [A"|) . 
and demonstrate explicitly that it reproduces the results 
expected in traditional MC simulations of these models 
f Section IIII A[) . In Section fill B[ we simulate extensions 
of the models, perturbed by weakening one bond per 
plaquette (an interaction that can arise experimentally 
in frustrated magnets when uniaxial pressure is applied 
}14{). In this case, the GLM update vastly outperforms 
traditional single-spin flips in the exploration of the low 
temperature physics of the models, in particular in find- 
ing a unique ordered groundstate. As discussed in Sec- 
tion IIV1 GLM update will help open the way for the ef- 
ficient simulation of a wide range of frustrated magnetic 
models in the future. 



II. GENERALIZED LOOP ALGORITHM FOR 
ISING MODELS. 

A. Ising models 

We consider classical Ising models that have Hamilto- 
nians of the general form 
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where spin variables Sf can take the value of ±1/2. In 
this paper, we refer to unperturbed models as those with 
Jij equal to a constant J > 0, while perturbed models 
(discussed in Section llll Bp have some of the Jij = J' / J 
(but still positive). The bond variables Sij are defined 
to have a value of either +1 ("AFM") or -1 ("FM") 
for each bond on the lattice. In a given spin configu- 
ration, a bond is referred to as satisfied if dijSfSj < 
and unsatisfied otherwise. An Ising model is referred to 
as frustrated if it is impossible for any spin configuration 
to satisfy all of the bonds on the lattice. We consider 
both perturbed and unperturbed versions of the follow- 
ing two-dimensional (2D) periodic lattice models in this 
paper (Fig [TJ : 

1. The triangular lattice AFM, where all Sij = 1. 

2. The fully frustrated square lattice. 

3. The fully frustrated honeycomb lattice. 

For the latter two lattices, uniform AFM interactions are 
unfrustrated; full frustration can be induced by enforc- 
ing the constraint that the product of the sign of the 
bond variables around each closed plaquette (square or 
hexagon) satisfies 
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-1, 



e.g. 



■1 for one bond, and <5, 



(2) 



1 for the rest. 



called fully-frustrated (FF). Thus defined, the three frus- 
trated Ising models that we are considering are known 
to have disordered groundstates with an extensive num- 
ber of equal-energy configurations jy, 0, Ho - Hq . The ge- 
ometry of the FM and AFM bonds for each lattice is 
illustrated in Fig. [T] 



B. Single-spin flips 

It is common to study the finite temperature (T) 
and groundstate properties of such models using Markov 
Chain Metropolis Monte Carlo methods. Traditionally, 
single spin flip (SSF) updates are employed in the Monte 
Carlo algorithm; one simply attempts to flip a single spin 
(one at a time) by computing the corresponding change 
in energy (AE) and then accepting the update with a 
Metropolis condition using probability: 



PfUp = min(l, e T ) 



(3) 



Since each plaquette is frustrated such interactions are 



In frustrated Ising models, many spin configurations are 
local energy minima, and most SSF attempts cost a large 
energy proportional to J. Based on Eq. ((3j), the proba- 
bility of accepting such a SSF update decreases expo- 
nentially as the temperature decreases. Consequently, at 
T <C J, once a SSF simulation has reached one of the lo- 
cal energy minima it is very unlikely for the simulation to 
accept any SSF updates. As a result, SSF updates tends 
to become "frozen" into one of the local energy minima 
at T -c J (i.e. ergodicity is lost). 

One can observe that in these local energy minima, 
there are often groups of spins that could be flipped to- 
gether without changing the number of unsatisfied bonds 
(and hence without significantly changing the energy of 
the system). Cluster algorithms take advantage of this 
observation to find and flip such groups of spins [8[ . 



C. A Generalized Loop Move 

Motivated by the above, we introduce a Generalized 
Loop Move Algorithm (GLM), designed to improve the 
sampling of equal (or nearly-equal) energy configura- 
tions near the groundstates of 2D periodic frustrated 
Ising models with bond sharing plaquettes. Like previ- 
ous al gori thms designed to work on specific Ising models 
B El, (HI j the algorithm works by identifying groups or 
clusters of spins that could be flipped together without 
changing the number of unsatisfied bonds, and then at- 
tempting to flip them according to a standard Metropolis 
update. In this way, no energy costs (proportional to J) 
are accrued, which as mentioned above can inhibit sim- 
ulation dynamics at low temperatures. 

To find these groups of spins, the GLM algorithm uses 
the concepts of dual graphs and matchings from graph 
theory. In this paper, dual graphs are defined in the 
usual graph theory sense. The dual nodes represent the 
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FIG. 1: Frustrated Ising spin models: AFM triangle (top 
left), FF square (top right) and FF honeycomb (bottom). 
The single bonds represent AFM bonds, the double bonds 
represent FM bonds, and the perturbed (weakened) bonds in 
the perturbed Ising models are indicated by "P" (see Section 

hitbT) . 



faces/plaqucttcs in the original ("primal") spin lattice, 
and the dual edges connect dual nodes which correspond 
to faces that share a bond in the primal spin system. See 
Fig. [2] for how the dual graphs are defined on the triangle, 
square and honeycomb lattices respectively. If we take a 
spin system and examine its dual graph, we can consider 
each dual edge to be inside a pseudo-matching when- 
ever the corresponding bond in the primal spin lattice 
is unsatisfied. This pseudo-matching is not a matching 
in the rigorous graph theory sense. In a graph theory 
matching, each node can only be adjacent to at most one 
edge in the matching. On the other hand, in the pseudo- 
matching described in this paper, each dual node can be 
adjacent to multiple dual edges in the pseudo-matching. 
There are a few interesting properties to note about this 
pseudo-matching: 

1. In any pseudo-matching that corresponds to a valid 
spin configuration, each dual node is adjacent to 
at least one dual edge in the pseudo-matching be- 
cause each plaquette contains at least one unsatis- 
fied bond in a fully frustrated Ising model. 

2. Each plaquette in a fully frustrated Ising model 
always has an odd number of unsatisfied bonds. 
Consequently in any pseudo-matching that corre- 
spond to a valid spin configuration, each dual node 
is adjacent to an odd number of dual edges in the 
pseudo-matching. 

3. Each valid ground-state spin configuration and 



its spin inverse corresponds to a unique pseudo- 
matching. 

4. There is a bijective correspondence between pairs 
of valid spin configurations and their spin inverses 
with the minimal number of unsatisfied bonds (one 
per plaquette) and perfect matchings in the dual 
graph. A perfect matching is a matching where 
each node is adjacent to exactly one edge in the 
matching. 

Given a spin configuration, the problem of finding an- 
other spin configuration with the same number of unsat- 
isfied bonds is equivalent to finding a pseudo-matching in 
the dual graph with the same number of dual edges as the 
pseudo-matching corresponding to the original spin con- 
figuration. In graph theory, a basic strategy for finding a 
new matching, with the same number of edges as an exist- 
ing matching, is to find an alternating cycle. Given some 
matching, an alternating cycle consists of edges that are 
alternatingly inside and outside the matching. For each 
edge in the alternating cycle, we can add it to the match- 
ing if it was not originally in the matching, and remove it 
from the matching otherwise, to get another matching of 
with the same number of edges. This operation is known 
as taking a symmetric difference between the matching 
and the alternating cycle. The same strategy can be 
used to find a new pseudo-matching with the same num- 
ber of dual edges as an existing pseudo-matching. Thus, 
if we start with a spin configuration and want to find 
another with the same number of unsatisfied edges, we 
can consider its corresponding pseudo-matching, use an 
alternating cycle to find another pseudo-matching, and 
then map this back to a spin configuration. If we look 
at this series of actions from the perspective of the spin 
system, we have found a group of spins delineated by the 
alternating cycle in the dual graph that can be flipped 
without changing the number of unsatisfied bonds. 

The motivation for looking at the dual graph of the 
spin system comes from the work by Bieche et. al. [l9(. 
In that work, the authors devised an algorithm that uses 
a modified dual graph of a spin system to find its ground- 
states. The algorithm takes the dual graph of the spin 
system and modifies it so that there is a bijective cor- 
respondence between perfect matchings in the modified 
dual graph and valid spin configurations. Weights are 
applied to each edge in the modified dual graph so that 
the weight of a perfect matching is equal to the energy 
of the corresponding spin configuration. Consequently, 
the problem of finding a groundstatc in the spin system 
is transformed into the problem of finding a minimally 
weighted perfect matching in the dual graph, for which a 
polynomial time algorithm has been devised by Edmonds 
et. al. |20,|2l[. One anticipates that the GLM update, 
designed in this paper to operate on fully frustrated Ising 
models, could be extended to work on a broader class of 
frustrated Ising models by adapting it to use the modified 
dual graph described in Biechc's paper. 




FIG. 2: Ising spin models and the corresponding dual graphs: 
AFM triangle (top left), FF square (top right), and FF hon- 
eycomb (bottom). The solid black lines represent primal spin 
lattice and the thick dashed red lines represent the dual lat- 
tice. The dashed shading of primal lattice along the verti- 
cal axis represent the Vertical Reference Cut and the hashed 
shading of primal lattice along the horizontal axis represent 
the Horizontal Reference Cut. 



D. Description of the Algorithm 

Like previous loop and cluster algorithms, The Gen- 
eralized Loop Move (GLM) is designed to supplement 
single spin flip (SSF) updates. Each iteration of the 
whole algorithm involves first attempting single spin flips 
on each spin in the lattice, and then attempting a fixed 
number of GLMs on the whole lattice. In this paper, we 
have chosen to perform N/30 GLM attempts after each 
SSF sweep (where N is the number of spins in the lat- 
tice). In each GLM attempt, we try to identify a group 
of spins that can be nipped together without changing 
the number of unsatisfied bonds in the lattice, compute 
the probability with which we should accept the flip so 
as to maintain detailed balance, and then accept the flip 
according to this probability. 

For ease of description, we divide the GLM algorithm 
into 3 subtasks: the Partition Finding subtask, the Loop 
Generating subtask and the Acceptance Probability Cal- 
culating subtask. The Partition Finding subtask de- 
scribes how the algorithm finds a group of spins that can 
be flipped together without changing the- number of un- 
satisfied bonds. The Loop Generating subtask describes 
how the algorithm finds cycles in the dual graph that 
are used by the Partition Finding subtask to identify the 
above-mentioned group of spins. Lastly, the Acceptance 
Probability Calculating subtask, as the name implies, de- 



scribes how the algorithm calculates the acceptance prob- 
ability with which we should accept suggested flips to 
maintain detailed balance. 

To describe the algorithm precisely, we will first estab- 
lish the following definitions and notations: 

1. A bond between spin i and spin j is referred to as 
Bij and its value is denoted as JijSij where J^ and 
Sij are defined as in Eq. ([!}; its dual is referred to 
as tijA . 



2. A set of dual edges is called a partition boundary 
if deleting all the bonds dual to these dual edges 
would divide the original lattice into two connected 
components. 

3. A Horizontal Reference Cut refers to a specific set 
of bonds that wraps horizontally around the lat- 
tice; a Vertical Reference Cut similarly refers to a 
specific set of bonds that wraps vertically around 
the lattice. These reference cuts are used in the 
algorithm to determine whether a set of dual edges 
constitutes a partition boundary. The set of bonds 
chosen is arbitrary. See Fig. [5] for how the Hori- 
zontal and Vertical Reference Cuts are defined in 
this paper on the AFM triangle, FF square and 
FF honeycomb lattices respectively. Other choices 
of horizontal and vertical reference cuts that wrap 
horizontally or vertically around the lattice would 
work just as well. 

4. A cycle is called a simple cycle if it does not inter- 
sect itself. If a simple cycle of dual edges in the dual 
lattice is dual to an even number of bonds in the 
Horizontal Reference Cut and an even number of 
bonds in the Vertical Reference Cut, then we refer 
to the cycle of dual edges as a closed loop. See Fig. 
[3] for an example of a closed loop. It can be shown 
that such a cycle of dual edges forms a partition 
boundary. 

5. If a simple cycle of dual edges in the dual lattice is 
dual to an odd number of bonds in the Horizontal 
Reference Cut or an odd number of bonds in the 
Vertical Reference Cut, then we refer to the cycle 
of dual edges as a cut. See Fig. [3] for an example 
of a cut. It can be shown that such a cycle of dual 
edges does not form a partition boundary. 

6. If the union of two non-intersecting cuts is dual to 
an even number of bonds in the Horizontal Ref- 
erence Cut and an even number of bonds in the 
Vertical Reference Cut, then it is called a comple- 
mentary union of cuts. See Fig. [3] for an example 
of a complementary union of cuts. It can be shown 
that such a union forms a partition boundary. 





a cut. If it is a closed loop, it forms a partition boundary 
and we have found the desired partition. If it is a cut, 
then we use the Loop Generating subtask to find a second 
cycle in the dual lattice. If the second cycle is another 
cut that does not intersect the first one and the union of 
the two cycles form a complementary union of cuts, then 
the union of two cuts forms a partition boundary and 
we have found the desired partition. If the second cycle 
does not meet the above requirements, then we abort the 
Generalized Loop Move. The reason why we want to find 
complementary unions of cuts is that in perturbed mod- 
els at T <C J, numerical results have shown that GLM 
using closed loops alone has a high probability of becom- 
ing frozen in a near-groundstatc. On the other hand, 
GLM using a combination of closed loops and union of 
cuts can efficiently sample from the near-groundstatcs. 
It should be noted that when the algorithm finds a cut 
and then a closed loop, it is possible to discard the cut 
and keep the closed loop as a partition boundary. 



FIG. 3: Types of cycles in the dual graph: closed loop (top 
left), cut (top right), complementary union of cuts (bottom). 
The dashed shading on the primal lattice along the vertical 
axis represents the Vertical Reference Cut and the hashed 
shading of primal lattice along the horizontal axis represents 
the Horizontal Reference Gut. The dashed green line on the 
dual graph represents the set of dual edges making up the 
cycle. The hash shaded area represents the spins enclosed by 
the Closed Loop or Complementary Union of Cuts. 



E. Partition Finding Subtask 

In the Partition Finding subtask, we try to find a par- 
tition boundary that is cither a closed loop or a comple- 
mentary union of cuts. Once we find such a partition 
boundary, we can use the Acceptance Probability Calcu- 
lating subtask (Section III Gp to compute the probability 
with which we should flip all the spins in the smaller of 
the two partitions formed by the partition boundary so as 
to maintain detailed balance. Then we flip all the spins 
inside the smaller of the two partitions with the probabil- 
ity just computed. When all the spins in either partition 
are flipped, all the satisfied bonds dual to the partition 
boundary become unsatisfied, all the unsatisfied bonds 
dual to the partition boundary become satisfied and all 
other bonds remain unchanged. To ensure that the num- 
ber of unsatisfied bonds remains unchanged when we flip 
all the spins inside the partition we want to choose a 
partition boundary that is dual to an equal number of 
satisfied and unsatisfied bonds. One way to satisfy this 
requirement is to choose a partition boundary that is al- 
ternatingly dual to satisfied and unsatisfied bonds. To 
find such a partition boundary, we use the Loop Gener- 
ating subtask (Section III F[) to find a simple cycle in the 
dual lattice that is alternating dual to satisfied and un- 
satisfied bonds. This simple cycle may be a closed loop or 



F. Loop Generating Subtask 

The Loop Generating subtask is designed to find sim- 
ple cycles of dual edges that are alternatingly dual to 
satisfied and unsatisfied bonds for the Partition Finding 
subtask to use. The main idea is that we start at a ran- 
dom dual node in the dual lattice, choose in a random 
but weighted way one of the dual edges adjacent to the 
dual node, go to the other dual node adjacent to the dual 
edge, choose in a random but weighted way a dual edge 
adjacent to this dual node (we will choose among dual 
edges dual to satisfied bonds if the last dual edge is dual 
to an unsatisfied bond and vice versa) and continue this 
way until we revisit one of the dual nodes we have passed 
through earlier (i.e. until the string of dual edges ends in 
a simple cycle). This simple cycle is used by the Parti- 
tion Finding subtask as described in Section Hi El More 
precisely, the algorithm can be described as follows: 

1. Randomly choose a dual node. Call it dualnodel 
and the set of dual edges adjacent to it E. 

2. Randomly select one of the dual edges Bij in E. 
The selection should be completely random for 
unperturbed models or based on some weighting 
scheme that favors strong unsatisfied bonds and 
weak satisfied bonds for perturbed models. When 
the spins inside the partition are flipped, the bonds 
dual to the partition boundary arc made satisfied 
if they were previously unsatisfied and unsatisfied 
if they were previously satisfied. So favoring strong 
unsatisfied bonds and weak satisfied bonds for in- 
clusion in the cycle means favoring partitions that 
would decrease the energy of the system if flipped. 
This increases the likelihood of finding a parti- 
tion that is energetically favored to be flipped at 
low temperatures. Numerical results show that a 



weighting scheme is particularly important for per- 
turbed triangle and square lattices where a GLM 
algorithm which selects partition boundaries com- 
pletely at random tends to fail to find partitions 



J 



that are energetically favored to be flipped in a 
reasonable amount of time. The weighting scheme 
chosen by the authors of this paper is: 
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if Bij is unsatisfied, 
if Bij is satisfied, 



(4) 



where f3 = 1/T and a is an arbitrary parameter 
that affects how strongly the weighting scheme fa- 
vors weak satisfied bonds and strong unsatisfied 
bonds j22[. The authors chose a to be 5. Pre- 
liminary numerical results showed that an a of 1 
did not sufficiently favor partitions that were ener- 
getically favored to be flipped. As a result, GLM 
with an a of 1 showed slow equilibration of thermo- 
dynamic data near the groundstate. Consequently, 
an a of 5 was tried and the thermodynamic data 
showed fast equilibration. 

3. Add this dual edge to our chain of dual edges, go to 
the other dual node adjacent to this dual edge and 
repeat the following steps until we revisit a dual 
node we have already visited: 

(a) Call the set of dual edges adjacent to the cur- 
rent dual node E 

(b) Remove from E all satisfied dual edges if the 
last dual edge added to our chain is satisfied; 
otherwise, remove all unsatisfied dual edges. 

(c) Remove from E all dual edges that would re- 
sult in a simple cycle with an odd number of 
dual edges if selected [23| ■ 

(d) If E is empty, then abort the Generalized Loop 
Move algorithm and leave the spin configura- 
tion unchanged. 

(e) Randomly select a dual edge from E based on 
a weighting scheme that favors dual edges dual 
to strong unsatisfied bonds or weak satisfied 
bonds. The author of this paper chose to use 
the weighting scheme defined by Eq. (QJ. 

(f ) Add the dual edge to the chain and go to the 
dual node on the other end of the dual edge. 

4. Once we revisit a node, the simple cycle at the end 
of our chain of dual edges is the desired simple cycle 
that we will use in the Partition Finding subtask. 

The computation time required to find a simple cycle 
through the algorithm described above is an important 
factor in the efficiency of the GLM algorithm. Numer- 
ical results concerning the amount of computation time 
required are discussed in Section IIIII 



G. Acceptance Probability Calculating Subtask 

Once the Partition Finding subtask generates a parti- 
tion boundary in the dual graph, the Acceptance Proba- 
bility Calculating subtask computes an acceptance prob- 
ability with which the GLM algorithm should flip all the 
spins inside the smaller partition so as to maintain de- 
tailed balance. For ease of exposition, it is convenient to 
establish the following definitions and notations: 

• Define a function / on the set of dual edge chains 
that end in a simple cycle as follows: If C = 
[ei e2 ... e a e a +\ ... e„_i e n ] is a chain of dual edges 
that ends in a simple cycle L = [e a e a +i ... e„_i e n ], 
then /(C) = [ei e 2 ... e n e„_i ... e Q+ i e a ]. Essen- 
tially, / takes a chain of dual edges that end in 
a simple cycle and reverses the order of dual edges 
that are part of the simple cycle (see Fig. S|). Note: 
f(f(C)) = C. 
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f(E) = 




FIG. 4: The function / described in the text. 

Let P se i ec t(C!, S) be the probability of the Loop 
Generating subtask selecting the dual edge chain 
C under the spin configuration S. 



• Denote the current spin configuration as Si and the Then the acceptance probability is defined by the follow- 
spin configuration after flipping all the spins inside ing expressions: 
the smaller partition as S 2 . 

I 



• If the Partition Finding subtask found a partition boundary that is a closed loop L from a chain C, then: 

PacceAC, Si) = min (l, e -m { S 2} -E {Sl)] Pselec t (f(C) S 2 ) \ (g) 

\ -r select (<-', <Jl] / 

In the unperturbed case where E(Si) = E(S 2 ), the expression simplifies to: 

Pselect(f(C), S 2 ) 



Pacce P t(C, Si) = min 1, 



Pselectiy, Si) 



(6) 



• If the Partition Finding subtask found a partition boundary that is a complementary union of cuts Li U L 2 from 
chains Ciand C 2 , then: 

iW((Ci,Ca), Si) = mm (l, e -g[ g (5 a )- g (50] ^(/(^). ^^^(/(Ca) g 2 ) \ (?) 

V "select (vi) ^>l)^select(^2, <->lJ J 



In the unperturbed case where E(S\) = E(S 2 ), the expression simplifies to: 

Pselect(f(Cl), S 2 )Pselect(f(C 2 ), S 2 ) 



Paccept((Ei,E 2 ), Si) = min 1, 



Pselect(Cl, Sl)P se l e ct(C 2 , Si) 
I 



(8) 



The authors of this paper have chosen to compute 
P S elect(C,Si) and P se i e ct(f(C),S 2 ) after a suitable par- 
tition boundary has been found. 

A rigorous proof of detailed balance is included in Ap- 
pendix |XJ 



III. SIMULATION RESULTS 

In this section we examine some prototypical data from 
Monte Carlo simulations of three frustrated Ising models, 
using simulations with SSF and GLM updates. In section 
IIIIA1 thermodynamic data for unperturbed Ising models 
is compared between simulations employing only SSF, 
and simulations using both SSF and GLM updates. Some 
GLM loop properties are characterized. In section UlIBl 
we examine the performance of the GLM updates on the 
three perturbed Ising models, which are constructed to 
have a finite-T phase transition to a ground state with 
a unique spin configuration. We demonstrate that GLM 
updates helps the simulation find the correct groundstate 
in cases where SSF updates fail. 



A. Results on unperturbed models 

Metropolis Monte Carlo simulations are performed on 
the three aforementioned Ising models: AFM triangular, 
FF square and FF honeycomb (Fig. [T]). A typical run 



employs an annealing technique where the simulation is 
begun at a high temperature (T = 0.5), and T is gradu- 
ally lowered through small steps until the system settled 
in its disordered degenerate groundstate. At each tem- 
perature, a fixed number of iterations of cither SSF or a 
combination of SSF and GLM is applied and measure- 
ment estimators such as energy and magnetization are 
collected. Each annealing simulation is performed twice: 
once with SSF only, and once with a combination of SSF 
and GLM updates. 

The first numerical result we consider is the internal 
energy of the Ising models under an annealing simulation. 
As illustrated in Fig. [SJ the internal energy measured as 
a function of temperature appears identical between SSF 
and GLM simulations, to within statistical error. This 
lends practical proof to the fact that the GLM updates 
are well behaved (do not disrupt detailed balance in the 
Monte Carlo procedure), and find the same degenerate 
manifold of states as the SSF updates alone. In addition 
to annealing simulations, we also tested the performance 
of SSF and GLM in quenched simulations where the en- 
tire simulation (including the warm-up or equilibriation) 
is run at a single temperature. We ran these single simu- 
lations at T = 0.05, finding that for the case of the AFM 
triangular and FF square lattice Ising models, both SSF 
and GLM were able to find a groundstate configuration, 
the efficiency gain through GLM not being significant in 
this case. However, for the FF honeycomb lattice, SSF 
simulations tend to be noticeably slow in finding the dis- 
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FIG. 5: Energy per spin (E/N) for the AFM triangular lattice 
(top), FF square lattice (middle), and FF honeycomb lattice 
(bottom). The theoretical true groundstate energy per spin 
for each system is -0.25 



ordered groundstate, a situation remedied by employing 
GLM updates. This difficulty for SSF alone to simulate 
the low-temperature thermodynamic behavior of the FF 
honeycomb lattice was also noted in the work by An- 
drews et. al. |15( where a cluster update specific to the 
FF honeycomb was proposed, fn Andrew's work, it was 
shown that in an annealing simulation, SSF alone tends 
to become frozen in one of the groundstates at low tem- 
peratures which is evident by the inability of the average 
magnetization to converge to zero (Fig. 3(c) of Ref. [Fa]). 
On the other hand, both the algorithm proposed in that 
paper and GLM can find the correct average magnetiza- 
tion through annealing or quenched simulations. 

A meaningful characterization of the low-temperature 
dynamics of the GLM updates can be determined via 
the study of the acceptance rate of the algorithm. We 
define the acceptance rate as the percentage of GLM at- 
tempts that result in successful group spin flips. The 
results are shown in Fig. [6] For the AFM triangular 
lattice, the acceptance rate increases as temperature de- 
creases and appears independent of lattice size. For the 
FF square lattice, the acceptance rate is high across the 
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FIG. 6: Acceptance Rate of GLM on the AFM triangular 
lattice (top), FF square lattice (middle), and FF honeycomb 
lattice (bottom). 



whole temperature range and appears to approach 1 as 
lattice size increases. For the FF honeycomb lattice, the 
acceptance rate is high across the temperature range, in- 
creasing as the temperature decreases, and apparently 
increasing with lattice size. Most importantly, the fact 
that the acceptance rate is large (of order unity) in all 
cases, particularly for T/J < 0.1, suggests that the GLM 
updates are successful in exploring the manifold of de- 
generate configurations that contribute to the disordered 
groundstate. This is in direct contrast to an algorithm 
employing SSF updates only, which are known to have 
exponentially suppressed acceptance rates for T <C J. 

The average number of dual edges that the GLM algo- 
rithm encounters in each successful attempt of the GLM 
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FIG. 7: Loop Size for the AFM triangular lattice (top), FF 
square lattice (middle), and FF honeycomb lattice (bottom). 
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FIG. 8: Chain Size for the AFM triangular lattice (top), FF 
square lattice (middle), and FF honeycomb lattice (bottom). 



(which wc will call chain size), and the average number of 
dual edges that form the partition boundary (which we 
will call loop size), provide a measure of the amount of 
work that the GLM algorithm needs to perform in each 
iteration. How these two values change with temperature 
and system size is of interest in determining the efficiency 
of the algorithm. As illustrated in Fig. for unperturbed 
FF Ising models, the loop size is independent of lattice 
size, and increases slightly as temperature decreases. As 
illustrated in Fig. [HJ for unperturbed FF Ising models, 
the chain size also appears to increase as temperature 
decreases. Chain size tends to increase with lattice size, 
however the rate of increase decreases with lattice size, 
so that the chain size appears to approach a limit as lat- 



tice size increases. The increase in chain size with lattice 
size for smaller systems is likely due to a finite-size ef- 
fect, meaning that for large lattices both chain size and 
loop size will saturate. Therefore, the work performed 
in each GLM attempt is proportional to a fixed multi- 
ple of the work performed in each SSF on a single spin. 
This computational complexity of the GLM attempt is 
the motivation behind each GLM iteration consistin g of 
one SSF sweep followed by a TV/30 GLM attempts [24[ 
in this paper. In this way, the work performed in each 
GLM iteration is expected to be comparable to a fixed 
multiple of the work performed in a SSF sweep. More rig- 
orous study is required to determine whether the compu- 
tational complexity of each GLM iteration is a multiple 
of the work performed in a SSF sweep. 
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FIG. 9: Energy per Spin for the Perturbed AFM triangu- 
lar lattice (top), Perturbed FF square lattice (middle), and 
Perturbed FF honeycomb lattice (bottom). 




B. Results on perturbed models 

Wc now explore the effect of anisotropy in the Ising in- 
teraction on the performance of the MC algorithm, with 
SSF and GLM updates. In the following, we define per- 
turbed Ising models by weakening one bond per plaquette 
in Eq. (fTJ), where two different strengths of Jij are used: 
J and J'. The location of weakened bonds are defined 
via the geometry illustrated in Fig. [TJ . The one weakened 
bond per plaquette has the strength J' / J = 0.9. In the 
AFM triangular case, Sij = 1 for all bonds, and weakened 
bonds are chosen arbitrarily to occur along the rows of 
the lattice; in the two FF models, the weakened bond 
is chosen to correspond to the FM bond with 8^ = — 1 
for simplicity. In each case, the degenerate manifold of 
groundstates is lifted resulting in two unique ground- 
states (related by spin inversion symmetry) where the 
unsatisfied bonds are placed uniquely at the perturbed 
bonds. This perturbation presents a greater difficulty for 
SSF MC algorithms to find the true groundstate. As il- 
lustrated in Fig. [9l SSF updates are unable to find the 
true groundstate energy in any of the three perturbed 
Ising models through an annealing procedure. GLM up- 
dates on the other hand arc able to find the true ground- 



FIG. 10: Acceptance Rate of GLM on the Perturbed AFM 
triangular lattice (top), Perturbed FF square lattice (middle), 
and Perturbed FF honeycomb lattice (bottom). 



state through annealing (as well as in quenched simula- 
tions ran at a single T -C J). It is interesting to note that 
at low temperatures, the internal energy curve for GLM 
is noticeably smoother than the internal energy curve for 
SSF alone. This suggests that GLM is able to converge 
to the equilibrium faster than SSF alone. 

As can be seen in Figure [T0l for the perturbed AFM 
triangular lattice model, the acceptance rate of GLM 
updates increases as the temperature decreases between 
T = 0.5 and 0.1, and then sharply drops off to be- 
low T = 0.1. For the perturbed fully frustrated square 
and honeycomb lattices, the acceptance rate is moder- 
ately high and stable between T = 0.5 and T = 0.05 
and sharply drops to below T = 0.05. This abrupt 
drop to zero corresponds to the finite-T phase transition 
expected to occur in these models at a temperature pro- 
portional to the perturbed energy scale, which the GLM 
updates are able to find cleanly. 

In the case of perturbed Ising models, temperature has 
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observation, the GLM algorithm tends to find local par- 
tition boundaries that do not grow with lattice size at 
all temperatures. As illustrated in Figure [T3J chain size 
behaves similar to loop size with respect to lattice size 
and temperature. These results on loop size and chain 
size suggest that in the groundstate of some perturbed 
Ising models (such as AFM triangle and FF square) the 
amount of work the GLM algorithm needs to perform in 
each GLM attempt will grow linearly with size of the lat- 
tice at low temperatures. This suggests that the amount 
of work done in each GLM iteration (one SSF iteration 
followed by x attempts of the GLM where x = N/30) will 
grow quadratically with system size at low temperatures. 
This expensive algorithm cost can be avoided in practi- 
cal situations by recognizing it as a sign of the underlying 
long-range order of the system - extensive sampling by 
the GLM update is not necessary in such a case. Regard- 
less, further simulation results are required to determine 
the computational complexity of the GLM update in this 
interesting case. 
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FIG. 11: Loop Size for the Perturbed AFM triangular lattice 
(top), Perturbed FF square lattice (middle), and Perturbed 
FF honeycomb lattice (bottom). 



a significant effect on loop and chain size. As illustrated 
in Figure [TT1 for perturbed AFM triangle and FF Square 
lattices, at high temperatures {T = 0.5 to 0.1) loop size 
appears independent of lattices size except at low tem- 
peratures (below the phase transition) where loop size 
appears to grow linearly with the shortest side of the 
lattice. This is expected based on observing GLM up- 
dates in action. For perturbed AFM triangle and FF 
square lattices at low temperatures, the GLM algorithm 
tends to encounter spin configurations where the only 
possible partition boundaries are those that wrap around 
the lattice. For the perturbed FF honeycomb lattice, 
loop size appears independent of lattice size regardless 
of temperature. This is also expected because, based on 



IV. DISCUSSION 

In this paper, we have presented a generalized loop 
move (GLM) update for a class of frustrated Ising models 
on two-dimensional lattices which are composed of bond- 
sharing plaquettes. The algorithm is designed to allow 
efficient Metropolis Monte Carlo updates in a degenerate 
or quasi-degenerate set of spin configurations that con- 
tribute to an extensive manifold of disordered low-energy 
states. We have thoroughly described the algorithm with 
the goal of making it easy to implement in a wide variety 
of models, and have provided a rigorous proof of detailed 
balance in the general case. 

We implemented and tested the algorithm on 
three prototypical lattice models, the antiferromagnetic 
triangular-lattice Ising model, and the fully-frustrated 
square- and honeycomb-lattice Ising models. All three 
models admit a disordered, extensively-degenerate man- 
ifold of configurations in their unperturbed groundstate. 
The GLM update is demonstrated to complement tradi- 
tional single spin-flip updates in allowing the simulation 
to find all configurations that contribute to the degener- 
ate manifold at low temperatures. In contrast to single 
spin-flips, the GLM update is able to sample configura- 
tions in the degenerate subspace with a very high ac- 
ceptance rate, which is independent of temperature for 
T <C J ■ The GLM achieves this while maintaining a very 
efficient O(N) complexity for each iV-site lattice. 

Perturbing the models by weakening one bond per pla- 
quette (triangle, square or hexagon) causes them to se- 
lect a unique spin configuration in their groundstate. Al- 
though single spin-flip updates typically cannot find this 
unique state due to the large energy barriers associated 
with sampling the quasi-degenerate manifold of states, 
we have shown that the GLM updates are capable of 
finding the equilibrium groundstate in a highly efficient 
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We expect that our current work in describing and 
characterizing the GLM updates will lead to their use 
in other frustrated Ising spin models of physical impor- 
tance. The GLM update, as described in this paper, 
can be straight-forwardly generalized to work on other 
frustrated Ising models on two-dimensional lattices of 
bond-sharing plaqucttcs. This includes models where the 
exchange JV,- is randomly frustrated, such as Edwards- 
Anderson spin glass models [lj| , which may facilitate the 
exploration of spin glass phases and phase transitions in 
a plethora of related models in the future. Finally, we 
expect that the GLM moves will be particularly useful in 
more realistic models of frustrated materials, where ad- 
ditional long-range interactions, local perturbations, or 
pressure-induced interaction anisotropy leads to the se- 
lection of a unique groundstate. 



FIG. 12: Chain Size for the Perturbed AFM triangular lattice 
(top), Perturbed FF square lattice (middle), and Perturbed 
FF honeycomb lattice (bottom). 
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Appendix A: Proof of Detailed Balance 

To prove detailed balance, we need to show that given 
any two spin configurations Si and 52, the following 
equation is satisfied: 



n(S!)r(S! -► s 2 ) = n(s 2 )T(s 2 -► Si) 



(Al) 



where II(S') denotes the Boltzmann probability of spin 
configuration S, and T(Si — > Sj) denote the transition 
probability in the Generalized Loop Move algorithm from 
spin configuration Si to spin configuration Sj. 



1. Preliminary notations and definitions: 

For ease of exposition, it is convenient to establish the 
following notations and definitions: 

• Let r be the set of all possible dual edge chains 
that the Loop Generating subroutine can generate 
(i.e. r is the set of dual edge chains that are al- 
ternatingly dual to satisfied and unsatisfied bonds 
and end in a simple cycle). 

• Take any 2 spin configurations Si and S 2 , and let: 

— Fi = {7 G F| 7 ends in a closed loop and flip- 
ping all the spins inside the closed loop takes 
the spin system from Si to 52} 

- fii = {(71,72) S (L x r)|7i,72 ends in cuts 
L\, L 2 and flipping all the spins inside L\\JL 2 
takes the spin system from S\ to S 2 } 



J 



Then, the proof for detailed balance follows naturally: 



- *i = riufii 



- T2 = {7 G r| 7 end in a closed loop and flip- 
ping all the spins inside the closed loop takes 
the spin system from S 2 to Si} 

- Q 2 = {(71,72) G (r x r)|7i,72 ends in cuts 
L\, L 2 and flipping all the spins inside Li UL 2 
takes the spin system from 62 to Si} 

- * 2 = r 2 u n 2 

- *i = {ip e *i| Paccepti^P, Si) > 0} where 
Paccept is defined as in Equation [5] or [7] 

- <F 2 = {ip G * 2 | Paccept(ip, S 2 ) > 0} where 
Paccept is defined as in Equation [5] or [7] 

- ^1 and ^2 are defined as above to facilitate 
the proof of detailed balance. 

- define F: T U (r x V) -> F U (T x T) such 
that F(r/) = /( 7 ) V 7 G r and F(( 7l , 72 )) = 

(/(71), /(72))V( 7 i,72)G(rxr) 

- Let F be the restriction of F to 1 i< 1 



For convenience, we will 

define P S e/ect((7i:72) ) S) 

Psdccthi, s)p select { l2 , s) v( 7 i, 72 ) g (r x r) 



also 



Note: 



T(Si 



S, 



H^m'. p seiect{^, Si) P acce pt (ip, Si), where 
(t, 3) = (1, 2) or (2, 1) 



2. Proof of Detailed Balance 

If we can prove the following 3 subclaims: 

1. i^('J'i) C >F 2 (this is used to prove subclaim 2). 

2. F is bijectivc between 1 3/ 1 and <F 2 . 

3. W> G *;, Ii{Si)P se lect{lp, Si) Paccept^, Si) 

U{S 2 )P se i ec t(F(^), S 2 )P accept (F(4,), S 2 ). 



S^ e ^' i II(S , i)P se i ec t(^, Si) Paccept^, Si) = S (/ , e ^' i n(S , 2)P se i ec t(F(V'), S 2 )P accept {F(lp), S 2 ) 

by subclaim 3 

^ e9 'Ji(Sl)Pselect(lp, Si) P acce pt (ip, Si) = Y, Se ^U(S 2 )P S elect(S, S 2 )P aC cept{S, S 2 ) 

because Fis bijective between "I^and \& 2 

^(Si)E^^> P se lect(tp, Si)P a ccept{4>, Si) = U(S 2 )T, S£ ^ P S elect(S, S 2 )P aC cept(S, S 2 ) 

U(Si)T(Si -> 5 2 ) = n(5 2 )T(5 2 -> Si) 
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3. Proof of subclaims 

Subclaim 1: F(#i) C ^' 2 

Take any ip e *i, show that F(ip) E * 2 
1. F(V>) e *2 because: 

(a) tf> £ V^i implies "0 ends in a closed loop or union of complementary cuts that forms a partition boundary 
around a group of spins that if flipped will take the spin configuration from Si to 5*2. Therefore, F(ip) ends 
in a closed loop or union of complementary cuts that forms a partition boundary around the same group 
of spins, which if flipped will take the spin configuration from 5*2 back to Si 

(b) and P se iect(F(ip), S 2 ) > because 

1p€ E *i => Paccepti^, Si)>0 

=► mm(l, e -ems 2 )-E {Sl) ] P'eiectm), 5 2 ) > 

PselectiV, Si) 



=> P select (F(i>), S 2 )>0 
2. 



PacceptiHlP), S 2 ) = mm(l, e -/8[g(g»)-g(g0] P select{?P, S 1 ) 

P{ W ' ' ^ PselectiFty), S 2 y 

because ^eticf^ Pseiecti^, Si) > 0. 
3. 1 & 2 =S> F(V>) £ * 2 - Therefore, #(*i) C * 2 

Subclaim 2: F is bijective between 9 1 and ^/ 2 

1. F is injective because i^(^x) — ^2 by subclaim 1 and take any ipi, ip 2 € ^i, F(ipi) = F(ip 2 ) =>• t/>i = ^ by 
definition of /() 

2. F is subjective because for any 6 G \]/ 2 , F(<5) G ^ by subclaim 1, and F(F(S)) = S. 

Subclaim 3: For any ip G *' 1; show that Tl(Si)Pseiect(ip, Si)P a cce P t{ip, Si) = n(S2)Pseiect(F(i/j), S2)Pacce P t(F(ip), S 2 ) 
Let a = Il(Si)P se i ec t(i), Si), and b = Tl(S 2 )P select (F(iP), S 2 ). Then, 

n(Si)P select (^, Si)P accept (<<p, Si) = IL(Si)P 3elect (iP, Si)mm(l, "^ 2 j ^^^ ) 

II(5l) Pselect{ip, Si) 

= a min(l, — ) 
a 

= min(a, b) 
U(S 2 )P select (FW, S 2 )P accept (F(V0, S 2 ) = U(S 2 )P select (FW, S 2 )min(l "^ /'^if'^] , ) 

11(^2] P se lect{F{ip), S 2 ) 

= b min(l, — ) 
= min(a, b) 

Therefore H(Si)P select {^, Si)P accept ^, Si) = Ii{S 2 )P select { F{jj), S 2 )P accept {F(4,), S 2 ) 



